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SECONDARY  INSTABILITY  IN  3D  MAGNETIC  RECONNECTION 


I.  INTRODUCTION 

Magnetic  reconnection  is  believed  to  be  responsible  for  a  wide  variety  of  energetic 
phenomena  in  both  laboratory  and  astrophysical  plasmas,  e.g.,  tokamak  disruptions,  ter¬ 
restrial  substorms,  and  solar  and  stellar  flares.  However,  in  astrophysical  situations  the 
magnetic  Reynolds  number  is  often  so  large  that  the  well-accepted  linear  and  non-linear 
reconnection  processes  are  too  slow  to  explain  the  phenomena  of  interest.  For  example,  it  is 
well  known  that  the  tearing  mode  growth  time  and  the  Sweet-Parker  reconnection  rate  are 
too  slow  to  account  for  solar  flares  1 .  Hence,  a  crucial  issue  for  astrophysical  plasmas  is  the 
determination  of  mechanisms  that  can  increase  the  rate  of  tearing.  Turbulence  is  probably 
the  most  commonly  invoked  process  by  which  energy  release  in  magnetic  reconnection  can 
be  increased  2  3  4  5  8. 

However,  the  way  in  which  a  laminar  reconnection  process  becomes  turbulent  has 
remained  somewhat  mysterious.  This  is  due  to  several  factors:  Simultaneous  observations 
of  all  the  relevant  fields  at  all  the  relevant  spatial  scales  are  presently  not  available  for  solar, 
space,  and  astrophysical  plasmas.  Laboratory  reconnection  experiments  are  very  difficult 
to  perform  and  diagnose,  and  theory  and  computation  are  hampered  by  the  complexity  of 
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the  governing  equations.  These  difficulties  are  not  as  severe  for  non-conducting  fluids,  and 
consequently  somewhat  more  progress  has  been  made  on  the  transition  problem  in  Navier- 
Stokes  fluids.  One  very  interesting  finding  has  been  that  for  a  wide  variety  of  simple  flows 
an  ideal  secondary  instability  can  be  the  first  step  in  the  transition  process  7  8.  When  the 
secondary  instability  attains  a  large  amplitude,  there  is  a  breakdown  of  the  laminar  flow 
structure,  generation  of  small  scale  structure,  and  chaotic  behavior  of  the  flow  variables. 

Paraphrasing  Orszag  and  Patera  9  the  three  steps  in  the  secondary  instability  process 

are: 

[1]  a  one- dimensional  primary  equilibrium  is  destabilized  by  a  two-dimensional  (2D)  pri¬ 
mary  linear  disturbance; 

[2]  the  primary  linear  disturbance  saturates  and  a  two-dimensional  secondary  equilibrium 
state  develops;  and 

[3]  the  two-dimensional  secondary  equilibrium  state  is  destabilized  by  an  ideal  three- 
dimensional  (3D)  secondary  instability. 

Important  flows  in  which  such  secondary  instabilities  occur  are  flat  plate  boundary  layers  9, 
plane  channel  flows  (Poiseuille  and  Couette),  10  9,  Hagen-Poiseuille  flow  9 ,  Taylor-Couette 
flow  11  and  free-shear  layers  12  1S.  It  is  well-known  that  the  tearing  instability  in  a  current- 
sheet  geometry  exhibits  both  step  [1]  and  step  [2].  Step  [1]  is  probably  the  most  extensively 
studied  problem  in  plasma  physics,  starting  with  the  basic  paper  by  Furth,  Killeen  and 
Rosenbluth  14 .  Step  [2]  was  first  demonstrated  by  Rutherford,  15  and  has  been  confirmed 
by  others  16 .  One  of  the  key  results  of  the  present  paper  is  that  reconnection  in  a  magnetic 
neutral  sheet  also  exhibits  Step  [3].  To  our  knowledge,  the  possibility  of  this  happening 
was  first  raised  by  Montgomery  17 . 

Why  is  a  three-step  process  necessary  for  increasing  the  rate  at  which  magnetic  energy 
is  released  ?  To  put  it  another  way,  since  going  to  three  spatial  dimensions  seems  to  make 
a  difference,  why  not  just  start  with  a  fully  three-dimensional  primary  perturbation  ?  This 
is  inadequate  for  the  following  reason.  The  two-dimensional  primary  linear  disturbances 
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grow  at  rates  which  are  too  slow  to  explain  the  phenomena  of  interest.  For  example, 
typical  growth  times  in  the  solar  corona  are  estimated  to  be  of  order  tens  of  days,  whereas 
solar  flares  are  observed  to  occur  on  the  time  scale  of  minutes.  In  section  2  of  this  paper, 
we  will  show  numerically  that  the  effect  of  increasing  three-dimensionality  of  the  primary 
modes  is  to  decrease  their  growth  rate.  This  also  can  be  shown  analytically  by  proving 
a  theorem  analogous  to  the  Squire’s  theorem  in  hydrodynamics,  which  states  that  under 
certain  very  broad  conditions,  two-dimensional  primary  perturbations  will  grow  faster 
than  three-dimensional  primary  perturbations  18 .  Hence,  the  time-scale  problem  exists 
a  fortiori  in  three-dimensions,  and  some  other  process  must  be  invoked  to  account  for 
fast  release  of  magnetic  energy.  In  this  paper  we  will  show  that  what  really  makes  a 
difference  in  timescales  is  the  basic  state,  i.e.  if  instead  of  the  standard  ID  neutral  sheet 
we  use  a  neutral  sheet  plus  a  large  dose  of  its  2D  unstable  primary  eigenfunction,  then 
ideal  instabilities  can  be  found. 

The  program  of  research  described  in  this  paper  requires  a  reformulation  of  some  of 
the  usual  MHD  theory  along  more  classical  hydrodynamic  lines,  elements  of  which  we  have 
reported  on  elsewhere.  The  problems  of  stability  and  transition  in  electric  current  sheets 
are  formally  similar  to  problems  in  hydrodynamics  19  17 .  For  example,  the  hydrodynamic 
free-shear  layer  is,  broadly  speaking,  a  configuration  similar  to  the  magnetic  neutral  layer 
13 .  This  suggests  that  techniques  developed  for  studying  the  hydrodynamic  transition 
problem  could  profitably  be  employed  in  investigations  of  magnetic  reconnection.  This 
has  been  done,  for  example,  for  primary  instabilities  of  magnetized  flows  20  21  22 . 

We  have  previously  derived  and  solved  a  2D  magnetohydrodynamic  analogue  of  the 
Orr-Sommerfeld  equation  23  for  an  investigation  of  Step  [1]  24 .  Two  significant  results 
emerged.  First,  it  appears  that,  for  linear  resistive  instability  to  occur,  it  is  necessary  that 
the  electric  current  profile  have  inflection  points.  Steep  gradients  are  not  enough.  Second, 
we  determined  the  existence  of  a  stability  boundary  which  depends  on  the  geometrical 
mean  of  the  viscous  and  resistive  Lundquist  numbers.  This  has  been  seen  by  others  25 
28 .  In  a  subsequent  investigation  of  Step  [2],  we  considered  the  nonlinear  evolution  of  two- 
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dimensional  perturbed  electric  current  sheets  27 .  We  observed  the  nonlinear  saturation 
of  unstable  eigenmodes.  The  secondary  equilibrium  state  was  close  to  that  predicted  by 
means  of  a  Landau  stability  equation.  This  analytical  technique  was  first  used  by  Stuart 
28  in  a  study  of  the  nonlinear  stability  of  plane  Poiseuille  flow  and  Taylor- Couette  flow. 
When  random  initial  conditions  were  used,  a  more  complex  ev  olution  toward  saturation 
occurred,  with  secondary  tearing  and  coalescence  observed  29  so. 

The  present  calculations  extend  the  previous  work  to  three  spatial  dimensions.  These 
calculations  differ  from  the  previous  ones  in  the  boundary  conditions  imposed  in  the  direc¬ 
tion  across  the  neutral  sheet.  Also,  they  differ  in  that  the  mean  magnetic  field  is  allowed 
to  diffuse,  i.e.  it  is  not  supported  by  an  electric  field.  Hence  the  saturated  state  of  the  2D 
disturbances  is  not  really  an  equilibrium  or  a  steady  state.  Instead,  the  diffusion  of  the 
mean  field  can  be  substantial.  However,  as  will  be  seen,  the  secondary  instability  discussed 
in  this  paper  seems  to  be  relatively  insensitive  to  this  decay. 

In  the  second  section  of  this  paper  we  set  up  the  foundation  for  a  discussion  of  the  fully 
three-dimensional  nonlinear  problem.  Here  we  first  describe  the  governing  equations,  and 
boundary  conditions.  The  primary  equilibrium  is  then  given,  followed  by  a  discussion  of 
its  two  and  three  dimensional  primary  resistive  instabilities.  An  Orr-Sommerfeld  equation 
is  derived  and  solved  for  the  primary  linear  instabilities.  Additional  Squire  equations  are 
given  for  obtaining  the  three-dimensional  primary  eigenfunctions.  These  previous  sections 
all  lead  to  the  third  section,  which  is  the  most  significant  part  of  this  paper.  In  this 
section  we  describe  the  evidence  for  secondary  three-dimensional  linear  instabilities.  We 
will  present  numerical  evidence  that  this  instability  is  ideal.  We  also  discuss  the  energetics 
of  the  mode  in  this  section.  In  particular,  we  will  show  that  the  dominant  energy  transfer  is 
from  the  one-dimensional  to  the  three  dimensional  fields.  We  will  further  show  that  onset 
of  the  instability  can  be  predicted  by  means  of  a  classical  potential  energy  analysis.  In  the 
fourth  section  we  give  a  preliminary  description  of  what  happens  when  the  secondary  modes 
attain  sufficient  amplitude  to  become  nonlinear,  using  high-resolution,  fully  nonlinear  3D 
numerical  simulations.  These  simulation  indicate  that  the  electric  current  sheet  becomes 
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turbulent,  and  then  “relaminarizes”  due  to  dissipation  into  a  new,  fully  three-dimensional 
tertiary  6tate.  We  will  report  on  these  simulations  more  fully  in  a  subsequent  paper.  The 
last  section  contains  some  concluding  remarks. 

II.  EQUILIBRIUM  AND  PRIMARY  INSTABILITY 
A.  The  Governing  Equations 

We  begin  with  the  dissipative,  incompressible,  magnetohydrodynamic  (MHD)  equa¬ 
tions,  written  in  a  dimensionless  rotation  form: 


&v  1 

—  =  v  x  u  -  VH  +  Ma  j  x  B  + 

(la) 

V  •  v  =  0, 

(16) 

dB  1 

=  VxvxB  -  — — V  x  V  x  B, 

&  Rm 

(1  c) 

V  •  B  =  0, 

(ID) 

where  v(x,<)  =  (u,v,w)  =  flow  velocity,  «(x,t)  =  V  x  v  =  vorticity,  B(x,  t)  =  magnetic 
field,  j(x,<)  =  V  x  B  =  electric  current  density,  H(x,t)  =  mechanical  pressure  +  kinetic 
energy  density,  R  =  Reynolds  number  (with  the  viscosity  assumed  to  be  constant  and 
uniform),  Rm  ~  magnetic  Reynolds  number  (with  the  resistivity  assumed  to  be  constant 
and  uniform),  and  M> i  =  Alfven  Mach  number  (  =  0  for  the  pure  Navier-Stokes  problem). 
In  this  representation,  B  is  measured  in  units  of  the  background  field  far  from  the  neutral 
sheet.  The  velocities  are  measured  in  units  of  the  Alfven  speed,  Ca  =  Bo/(4irp)i  where 
the  mass  density,  p,  is  assumed  to  be  constant  and  uniform.  Hence,  Ma  =  1.  For  a 
characteristic  distance,  Lo ,  defined  by  the  neutral  sheet  width,  the  time  is  measured  in 
units  of  the  Alfven  transit  time,  Lq/Ca •  The  regimes  of  applicability  of  such  incompressible 
models  have  been  discussed  elsewhere  31  32  33 . 
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B.  Boundary  Conditions 

In  analogy  with  the  hydrodynamic  problems,  we  refer  to  the  spatial  coordinate  along 
the  magnetic  field  (x)  as  the  fieldwise  direction,  the  spatial  coordinate  across  the  magnetic 
neutral  sheet  (2)  as  the  cross-sheet  direction,  and  the  remaining  spatial  coordinate  (y)  as 
the  sheetwise  direction.  We  consider  a  system  with  periodic  boundary  conditions  in  the 
fieldwise  and  sheetwise  directions.  In  the  cross-sheet  direction  we  want  the  magnetic  field 
to  merge  into  the  background  field  as  it  moves  away  in  2  from  the  neutral  sheet.  Hence 
we  want  to  enforce 

Bx{±  00)  =  ±1;  By(±oo)  =  Bz(±  00)  =  0.  (2) 

We  impose  this  computationally  with  free-slip  boundary  conditions  34 : 


dBx 

dz 


dBy 

dz 


=  Bz  =  0  at  Lz 


(3«) 


and 


du 

dz 


—  =  w  =  0  at  Ls, 
dz 


where  L3  is  the  computational  boundary  in  the  cross-sheet  direction. 


(36) 


C.  Primary  Equilibrium 


To  obtain  the  primary  instability,  we  assume  the  following  linearized  decomposition: 


v(x,y,z,f)  =  U0{z)ex  +  v'(z,y,z,f), 


(4a) 


B(x,y,2,f)  =  #0(2)6*  +  b'(x,y,z,f), 


(46) 


and  for  the  pressure: 


P(x,y,z,t)  =  P0{z)  +  p'(x,y,z,t). 


(4c) 
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For  the  primary  neutral  sheet  equilibrium  we  choose 


Uo(z)  =  0,  (5a) 

and, 

i?o(z)  =  tanh  z.  (56) 

The  pressure  is  determined  by  taking  the  divergence  of  equation  la,  and  solving  the  re¬ 
sulting  Poisson  equation.  This  configuration  has  been  studied  by,  among  others,  Schnack 
and  Killeen  35 . 

D.  Primary  Two-  and  Three-dimensional  Instabilities 

After  linearizing  equations  la  —  Id  about  the  equilibrium  given  in  equations  5,  we 
assume  that  the  first-order  terms  can  be  decomposed  in  the  following  manner: 

/'(*,y,z,t)  =  f(z)  c»a*+‘0v-***'<,  (6) 

in  which  a  is  the  fieldwise  wavenumber,  /?  is  the  sheetwise  wavenumber,  and  u>  is  the 
complex  growth  rate.  After  eliminating  the  pressure,  the  primary  linear  perturbations  we 
will  employ  are  unstable  eigenfunctions  of  the  following  equations,  where  D  =  d/dz  and 
the  tildes  have  been  dropped: 

{ D 2  -  (a2  +  02)}2w  -  iaRU0{D 2  -  (a2  +  02)}w  +  iaR(D2U0)w  = 

(7a) 

-  iu>R{D2  -  (a2  +  02)}w  +  iaRM2A[{D2BQ)bz  -  B0{D 2  -  (a2  +  02)}bz], 

{D2  -  (a2  +  /?2)  -  iaRmUo  +  iu>Rm}bz  =  -iaRmB0w,  (76) 

with  the  boundary  conditions: 

w  ±  (oo)  =  Dw(± oo)  =  6z(±oo)  =  0.  (7c) 
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Equations  7  are  a  generalization  of  the  3D  Orr-Sommerfeld  equation  to  include  magnetic 
effects.  Note  also  that  we  have  written  equation  7  in  a  general  form  which  includes  a  mean 
flow. 

The  linear  equations  are  solved  using  the  Chebyshev  r  technique,  which  was  first 
applied  to  the  Orr-Sommerfeld  equation  by  Orszag  36  The  mesh  is  stretched  in  the  ver¬ 
tical  direction  with  a  hyperbolic  mapping  13 .  For  the  velocity  free  shear  layer,  the  code 
reproduces  the  linear  growth  rates  reported  by  Metcalfe  et  al.  13  when  Ma  =  0. 

It  can  be  shown  that  in  the  limit  of  zero  viscosity  and  resistivity,  there  are  no  unstable 
modes  of  equation  7.  A  proof  of  this  point  is  given  in  Appendix  A.  When  the  resistivity 
and  viscosity  are  finite,  the  unstable  modes  of  equation  7  have  an  imaginary  eigenvalue, 
real  w,  and  imaginary  bz.  Growth  rates  for  some  of  these  modes  are  provided  in  Table  I, 
where  a  unit  magnetic  Prandtl  (i.e.  Rm  =  R)  number  has  been  assumed.  All  of  the  cases 
listed  in  Table  I  have  a  =  0.2.  Note  that  the  system  size  is  scaled  by  width  of  the  neutral 
sheet,  so  that  a  =  0.2  implies  a  wavelength  for  the  perturbation  of  107r  times  the  sheet 
width.  In  Table  I,  the  2D  primary  instability  growth  rates  are  listed  under  J2D,  and  the 
3D  rates  are  listed  under  73 p.  The  last  column  of  Table  I  will  be  discussed  in  the  section 
on  secondary  instability.  Of  particular  interest  is  the  well-known  decrease  in  the  growth 
rate  as  the  resistive  Reynolds  number  increases  (note  that  since  we  have  set  R  =  Rm,  the 
growth  rates  in  Table  I  also  reflect  viscous  effects,  i.e.,  a  stabilization  as  R  decreases). 

Figure  1  shows  some  eigenfunctions  for  this  problem  to  illustrate  the  effect  of  increasing 
the  sheetwise  wavenumber.  In  this  figure,  the  modes  are  normalized  to  the  maximum  value 
of  bz.  Figure  la  is  the  familiar  primary  2D  unstable  eigenfunction,  while  figures  lb  and 
lc  show  unstable  modes  with  0  values  of  0.2  and  0.5,  respectively.  The  mode  appears  to 
concentrate  about  the  neutral  line  with  increasing  /?,  and  the  ratio  of  the  maximum  of  w  to 
the  maximum  of  bz  decreases.  Figure  Id  shows  a  purely  damped  mode  for  (3  =  2.  We  note 
here  that  the  linear  analysis  also  can  be  done  for  more  physically  interesting,  i.e.  much 
higher,  values  of  the  Reynolds  numbers.  We  limit  ourselves  here  to  Reynolds  numbers 
for  which  we  can  presently  perform  accurate,  three-dimensional,  fully  nonlinear  numerical 
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simulations. 


The  most  important  result  for  this  paper  concerns  the  variation  of  primary  instability 
growth  rates  with  sheetwise  wavenumber.  Figure  2  shows  some  computed  unstable  eigen¬ 
values  of  equations  7  as  a  function  of  the  sheetwise  wavenumber,  0.  It  is  evident  from  the 
figure  that  increasing  0  leads  to  a  decrease  in  the  primary  instability  growth  rate.  Above 
/ 3  =  0.2,  there  is  a  linear  decrease  in  growth  rate.  At  high  enough  0  a  complete  stabiliza¬ 
tion  occurs.  This  can,  in  fact,  be  shown  analytically  by  proving  a  Squire’s  theorem  for  the 
3-D  system,  *.e.,  the  three  dimensional  linear  modes  always  have  a  lower  growth  rate  than 
the  2D  modes.  For  this  reason,  one  might  conclude  that  3D  modes  could  be  safely  ignored 
in  considering  the  evolution  of  this  system.  The  calculations  presented  in  section  3  will 
invalidate  this  assumption. 

Once  the  eigenfunctions  have  been  determined  from  equations  7  the  other  components 
of  the  fields  must  be  determined  for  use  in  initializing  the  nonlinear  code.  In  two  spatial 
dimensions  this  is  easily  done  by  means  of  the  solenoidality  relation  for  the  velocity  and 
magnetic  fields.  In  three  spatial  dimensions,  we  must  solve  a  perturbed  vorticity  (or  Squire) 
equation: 


{ D 2  —  (a2  +  02)  —  iaRUo}(0u  —  av)  +  iaM2A  RBo(0bx  —  aby)  = 

—  iu?R(0u  —  a«)  —  (3R[(DUo)w  —  M\  ( DBo)bz ]. 
and  a  similar  equation  for  the  perturbed  electric  current: 


(8a) 


{D2  -  (a2  -f  02)  -  iaRjnUo}(0bx  -  aby )  -f-  iaRmBo(0u  -  av)  = 

(8b) 

-  iu)Rm(@bt  -  aby)  -  0Rm[(DUo)bz  -  ( DB0)w ]. 
with  the  boundary  conditions: 

/?u(±oo)  -  qv(±oo)  =  0  (8c) 
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(3bx(±  oo )  —  a6v(±oo)  =  0 


m 


The  remaining  field  components  then  can  be  determined  by  using  the  solenoidality  relations 
for  the  velocity  and  magnetic  field. 

E.  Secondary  Equilibrium  -  Nonlinear  Saturation  of  2D  Primary  Instability 

Since  the  2D  unstable  eigenmode  has  the  fastest  growth  ^ate,  it  is  reasonable  to 
suppose  that  it  will  be  the  first  mode  to  develop  to  the  point  where  nonlinear  effects 
become  significant.  It  has  been  determined,  both  analytically  and  numerically,  that  this 
mode  saturates  into  a  secondary  equilibrium  when  it  becomes  sufficiently  large.  Rutherford 
15  in  his  classic  analysis  showed  that  the  tearing  mode  changes  from  exponential  to  linear 
growth  at  a  mode  amplitude  which  varies  approximately  inversely  with  magnetic  Reynolds 
number.  Hence,  the  mode  saturates  at  very  low  amplitude,  whereupon  the  reconnection 
rate  becomes  that  due  simply  to  ohmic  diffusion.  The  conclusion  is  that  tearing  in  two 
dimensions  cannot  account  for  the  rapid  energy  release  observed  in  phenomena  such  as 
solar  flares. 

Figure  3  shows  the  magnetic  topology,  electric  current,  and  vorticity  for  a  system  close 
to  the  saturated  state.  This  resembles  the  initial  state  that  we  use  for  the  three-dimensional 
perturbation.  We  note  in  Figure  3a  the  familiar  magnetic  island  \riiich  characterizes  the 
tearing  mode,  in  Figure  3b  the  electric  current  filamentation  at  the  magnetic  X-point  37 , 
and  in  Figure  3c  the  flow  pattern.  Typically,  the  kinetic  energy  of  the  2D  saturated  state 
is  much  less  than  the  magnetic  energy,  which  demonstrates  that  the  evolution  has  become 
dominated  by  diffusion. 

III.  SECONDARY  INSTABILITY 

In  this  chapter  we  provide  numerical  evidence  for  the  existence  of  a  secondary  insta¬ 
bility  of  the  neutral  sheet,  investigate  its  parametric  dependence,  and  show  that  the  onset 
of  the  instability  is  predictable  by  a  classical  potential  energy  analysis.  The  investigation 
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relies  on  numerical  simulation  of  Equations  1.  The  semi-implicit  Fourier  pseudospectra!  - 
Fourier  collocation  scheme  that  we  employ  is  described  in  Appendix  B. 

A.  Evidence  for  a  secondary  instability 

We  initialize  the  magnetic  field  in  the  following  way: 

B  (*,»,*,<  =  0)  =  £0(z)e,  +  e2dB2d(z)zia*  +  tidB  sd(z)eia^+i^ ,  (9a) 

where  t2d  is  the  amplitude  of  the  two-dimensional  MHD  Orr-Sommerfeld  eigenmode  (B2(f), 
e2d  is  the  amplitude  of  the  three-dimensional  MHD  Orr-Sommerfeld  eigenmode  (Bjd),  and 
<}>  is  the  phase  shift  between  the  2D  and  the  3D  mode  (set  to  7r/2  throughout).  For  the 
velocity  field  we  have: 

\{x,y,  z,t  =  0)  =  e2dv2d(z)eia *  +  tid\id{z)exa^*)^v .  (96) 

The  magnetic  field  eigenfunctions  are  normalized  so  that  m.axz\Bz(z)\  =  1,  and  the  velocity 
eigenfunctions  are  scaled  proportionately  (see  Figure  1). 

We  now  demonstrate  that  the  nonlinearly  saturated  two-dimensional  tearing  layer 
is  subject  to  three-dimensional,  secondary  instabilities.  To  do  this  we  first  prepare  an 
initial  state  consisting  of  the  tearing  layer  plus  a  large  dose  of  the  2D  primary  unstable 
eigenfunction  (e2<*  =  -02),  so  that  the  system  is  close  to  the  2D  nonlinearly  saturated  state. 
Previous  research  has  shown  that  the  form  of  the  perturbation  is  relatively  unchanged 
in  the  saturated  state,  i.e.,  the  “shape  assumption”  holds  27  28 .  We  then  perturb  this 
state  with  a  much  smaller  dose  of  the  appropriate  3D  primary  unstable  eigenfunction 
(e3<£  =  1  x  10-8).  These  initial  conditions  represent  a  state  that  is  highly  probable,  since 
the  primary  2D  modes  outstrip  the  3D  primary  modes  and  saturate  first,  creating  a  kind 
of  secondary  equilibrium  in  which  the  3D  disturbances  subsequently  evolve. 

One  respect  in  which  our  simulation  might  differ  from  others  is  that  we  do  not  hinder 
the  evolution  of  the  current  layer,  but  instead  allow  it  to  dissipate  freely.  We  can  quantify 
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the  spreading  of  the  current  layer  by  defining  a  mean  current  thickness  as: 

tc  =  Bext[- r<-^-)]moi>  (10) 

where  V  equals  the  system  volume,  and  Bext  denotes  the  magnetic  field  exterior  to  the 
tearing  layer,  *.e.,  1.  Furthermore,  the  bracket  notation  is  used  to  denote  volume  integrals: 
(/)  =  /*°o  /  dx  dy  dz.  Initially,  lc  =  1  by  definition.  For  case  9  (see  Table  I) 

with  Rm  =  200,  at  t  =  200,  ie  =  3.10.  At  the  same  time,  for  case  18  with  Rm  =  1000, 
te  —  1.48.  Hence  the  spreading  of  the  current  layer  due  to  dissipation  is  significant. 
However,  as  we  will  show,  the  evolution  of  the  secondary  mode  appears  to  be  relatively 
insensitive  to  this  spreading. 

Although  we  initialize  the  fields  with  one  mode  in  x,  after  a  few  Alfven  times  when 
the  system  is  near  the  saturated  state,  more  than  one  mode  will  be  excited.  Hence,  it  is 
more  accurate  to  represent  the  fields  as  having  the  form  of  a  2D  tearing  layer  perturbed 
by  a  3D  disturbance: 

b  =  b'»(z,()*.+  e  b 2(MK“m'+x;  e  (lla) 

|m|<Af/2  n=±l  |m|<M/2 

and 


|m|<M/2  n=±l  |m|<M/2 

where  the  tilde  denotes  the  Fourier  coefficient. 

To  track  the  evolution  of  the  3D  secondary  modes,  we  define  a  3D  total  energy  given 
by 

E,  D=\j  (|B(s>r  +  |v<s>|V*.  (12a) 
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This,  as  well  as  the  2D  energy,  given  by: 

=  (|B<2>P  +  |v(2)|2)rfsi,  (12  b) 

is  shown  in  figure  4,  for  case  9  of  Table  I.  After  about  fifty  Alfven  times  E^d  flattens  out 
and  a  linear  phase  is  evident  in  Esd-  As  previously  noted,  the  growth  of  the  3D  energy  is 
relatively  independent  of  the  dissipative  spreading  of  the  neutral  sheet.  Results  of  several 
of  these  runs,  with  different  values  of  /?  and  the  Reynolds  numbers,  are  shown  in  Table  I. 
In  this  table  we  give  the  total  3D  energy  average  growth  rates  for  the  secondary  instability, 
7,  where 


1  d,E$D 
2  E5D  dt  ~7, 

and  200  <  t  <  250.  For  case  9,  for  example,  the  secondary  mode  grows  at  more  than 
double  the  rate  of  the  primary  3D  mode.  All  of  the  runs  listed  in  Table  I  were  performed 
with  a  =  0.2.  Limited  tests  of  a  variation  for  case  3  in  the  table  indicated  a  falloff  in  7  for 
higher  and  lower  a  than  0.2.  We  did  not  investigate  further  this  rather  fortuitous  result. 

Several  conclusions  can  be  drawn  from  the  information  presented  in  Table  I:  (a.)  There 
is  an  increase  in  the  secondary  mode  growth  rate  at  low  Reynolds  number,  then  a  relative 
independence  at  higher  Reynolds  numbers;  (b.)  There  is  an  instability  threshold  in  the 
Reynolds  number,  i.e.,  the  secondary  mode  can  be  stabilized  if  there  is  enough  diffusion. 
This  might  explain  why  this  mode  is  not  observed  in  the  reported  collisionless  tearing 
experiments,  which  have  had  magnetic  Reynolds  numbers  of  10  to  20  38  39 ,  and  (c.)  There  is 
a  high  /?  cutoff  of  the  secondary  instability.  The  cutoff  value  in  /?  increases  as  the  Reynolds 
numbers  increase  and  is  given  approximately  by  (3C  ~  .25 Rm  ■  However,  this  estimate  is 
clouded  somewhat  by  our  use  of  a  unit  magnetic  Prandtl  number.  Furthermore,  it  is  of 
interest  that  linearly  stable  3D  primary  eigenmodes  can  trigger  3D  secondary  instabilities. 
Unfortunately,  we  are  somewhat  limited  in  the  magnitude  of  Reynolds  numbers  that  we 
can  simulate  accurately  with  direct  numerical  simulation.  However,  we  note  that  even  at 
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these  modest  values,  the  growth  rates  for  the  secondary  instability  can  exceed  the  resistive 
primary  instability  growth  rates.  In  addition,  these  latter  rates  should  decrease  as  the 
magnetic  Reynolds  number  is  increased. 

B.  Energetics  of  the  secondary  instability 

While  the  behavior  of  7  is  indicative  of  a  secondary  linear  instability,  it  provides  little 
insight  into  the  dynamics  of  the  unstable  mode.  In  particular,  we  would  like  more  infor¬ 
mation  on  the  interactions  between  the  ID,  2D,  and  3D  fields.  To  obtain  this  information 
we  generalize  the  hydrodynamic  analysis  of  Orszag  and  Patera  (1983)  to  the  MHD  case. 

Hence,  we  now  decompose  the  magnetic  and  velocity  fields  so  that  the  various  fields 
are  mutually  orthogonal  (time  dependence  of  the  fields  is  assumed):  The  one  dimensional 
fields  are  given  by: 


bW(2)  =  (£<*>(*)  +  ?(*))«.  , 


(14o) 


and 


vW(z)  =  (£7(1)(z)  +  vj?(*))e.  =0  +  0  =  0  (146). 

For  the  two-dimensional  fields  we  have: 

BN(x,  z)  =  B<!>(x,  *)  -  B[l\z)e, ,  (14c) 

vlJl(x,2)  =  v(2)(x,2)  -  t>oo+)®**  (14<i) 

The  three-dimensional  fields  are  given  by: 

B!s](x,y,z)  =  B(3)(x,y,z),  (14e) 

v[3](x,y,z)  =  v(3)(x,y,z).  (14/) 
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The  total  energies  (magnetic  plus  kinetic)  corresponding  to  these  fields  evolve  as: 


=  —  T12  —  Tu  +  D  i, 


(15a) 


=  Tu  —  Th  -I-  Di, 


(15*) 


=  Tn  +  Tn  +  Dj. 


(15c) 


Our  interest  is  in  the  last  equation,  where  Ei  represents  the  energy  in  the  3D  fields, 
Tn  represents  the  interaction  of  the  ID  and  3D  fields,  Tn  represents  the  interaction  of 
the  2D  and  3D  fields,  and  Di  represents  the  ohmic  and  viscous  diffusion  of  the  3D  field. 
Interpretation  of  the  other  terms  follows  from  this. 

It  is  straightforward  to  evaluate  Tu  from  the  equation  15a,  viz., 


Til  =  -  J  BM  •  (V  x  vM  x  B^)dsx. 


T23  is  then  given  by: 


Tu  =  J  vM  •  (v  x  u)d3x  +  J  vM  •  (V  x  B  x  B)d3x  +  J  B(s]  •  (v  x  B)rf3x  -  T 13,  (17) 


and  Di  is  given  by: 

Di  =  J  | V  x  BW|2dsx  +  i  J  |V  x  vf3)|2d3x.  (18) 

To  obtain  a  rate  for  energy  transfer  to  and  from  the  3D  mode,  we  compute 

+  k  +  k  =  +  <19) 

where^j  =  J'(|Bl*J|2  +  |vi3l  )2)dsx.  To  check  the  accuracy  of  this  decomposition  we  compare 
7  with  7. 
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Figure  5  shows  Tu,T2s,  Ds  and  7  as  functions  of  time  for  case  9  of  Table  I.  It  is 
apparent  that  the  total  growth  rate  is  almost  entirely  due  to  the  interaction  between  the 
ID  and  3D  fields.  The  2D  field  extracts  energy  from  the  3D  field  during  the  secondary  linear 
phase.  The  dissipation  of  course  provides  an  energy  sink  throughout  the  run.  Figure  6 
shows  the  same  quantities  for  case  11,  to  show  how  the  secondary  mode  transfer  terms  vary 
with  /?.  This  case  is  near  to  the  maximum  growth  rate  observed  in  /?,  all  other  parameters 
being  the  same  (see  Table  I).  An  increase  in  T\ j  is  seen,  but  this  is  almost  balanced  by 
a  decrease  in  Ti$.  Figure  7  shows  the  same  quantities  for  case  13,  which  is  nearly  stable. 
It  can  be  seen  that  increases  in  (3  lead  to  increased  dissipation  through  decreased  D 3, 
whereas  the  transfer  term  T js  does  not  change  very  much.  Finally,  Figure  8  shows  the 
same  quantities  for  case  16,  to  illustrate  the  effect  of  increasing  the  Reynolds  numbers. 
Comparing  this  to  figure  5,  we  see  that  T 13  remains  about  the  same,  Dj  increases,  and 
T 23  decreases.  However,  the  secondary  mode  growth  rate  remains  about  the  same. 

C.  Potential  Energy  Analysis 

As  previously  noted,  the  primary  equilibrium  is  stable  to  ideal  perturbations  (we 
prove  this  point  in  Appendix  A).  On  the  other  hand,  we  have  also  seen  evidence  that 
the  secondary  instability  is  ideal.  The  secondary  instability  growth  rates  listed  in  Table  I 
appear  to  be  relatively  insensitive  to  changes  in  the  Reynolds  numbers.  Detailed  analysis  of 
energy  transfer  shows  that  the  ID  to  3D  transfer  also  is  relatively  unchanged  by  variation  of 
Reynolds  numbers  (compare  figure  5  and  8).  These  observations  suggest  that  the  stability 
of  the  secondary  mode  might  be  predictable  by  means  of  a  classical  energy  analysis  for 
ideal  instability  40 . 

Since  the  kinetic  energy  of  the  secondary  equilibrium  is  much  less  than  the  magnetic 
energy,  we  assume  that  the  secondary  equilibrium  state  is  static,  i.e., 

v<»>  =  0,  (20a) 

and 

=  B^(z)ez  +  BW(x,  z).  (206) 
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We  can  then  calculate  the  potential  energy  8W  relevant  to  our  system 


SW  =  ^-  J{ |V  X  [vM  x  B<‘>]|J  -  [v|Ji  X  (V  X  B<‘>)]  ■  (V  x  (vW  x  B<‘>J)}<lV 


(21) 


Note  that  the  surface  terms  drop  out  due  to  our  boundary  conditions.  If  6W  is  less  than 
zero,  the  system  is  unstable.  Figure  10  shows  the  same  data  as  in  figure  4,  along  with 
the  sign  of  6W  as  a  function  of  time  for  case  9  of  Table  I.  Note  that  6W  remains  positive 
approximately  until  the  2D  mode  has  saturated.  It  then  becomes  negative  at  about  the 
time  that  the  3D  mode  takes  off.  This  provides  additional  evidence  that  the  secondary 
instability  is  ideal,  i.e.,  independent  of  the  Reynolds  numbers. 


IV.  EVOLUTION  OF  PERTURBATIONS  INTO  THE  FULLY  NONLINEAR 
REGIME 

One  important  remaining  question  is  what  happens  when  the  secondary  modes  become 
nonlinear  ?  Is  there  a  transition  to  turbulence,  or  does  the  mode  saturate  into  a  tertiary 
equilibrium  ?  Here  we  briefly  describe  some  high  resolution  runs  which  suggest  that  the 
answer  is  both,  i.e.  there  is  a  burst  of  turbulence  followed  by  a  relaminarization  due  to 
the  relatively  high  dissipation  employed  in  the  simulation.  The  runs  are  initialized  with 
eigenfunctions  which  evolve  into  the  nonlinear  regime,  with  a  =  0.2,  (3  =  0.2,  R  =  Rm  = 
50, 100,  and  200,  C2 d  =  0-01  and  =  0.001.  The  resolution  for  these  runs  is  (32  x  32  x  64) 
modes.  Some  data  for  these  runs  are  given  in  Table  II. 

Evidence  for  a  burst  of  turbulence  is  provided  by  an  examination  of  the  enstrophy 
for  cases  1-3  in  Table  II.  The  kinetic  energy  decay  rate  is  proportional  to  the  enstrophy, 

ft  =  <M2)>  *-c., 

^  =  -i(M2)  (22) 

and  it  can  be  used  as  a  proxy  for  monitoring  the  level  of  small  scale  excitation.  This 
quantity  is  shown  for  three  runs  at  different  Reynolds  numbers  in  figure  10.  For  all  of  the 
cases  there  is  a  big  jump  in  fl  around  t  =  100,  followed  by  a  slow  decrease.  The  magnitude 
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of  the  burst  increases  with  the  Reynolds  numbers,  as  does  the  subsequent  level.  Similar 
behavior  is  seen  in  the  perturbed  mean  square  electric  current. 

The  evolution  is  clarified  by  the  time  histories  of  the  perturbed  components  (note 
the  normalization).  Figure  11  and  12  show  the  logarithm  of  selected  modal  kinetic  and 
magnetic  energies,  with  the  z  direction  integrated  out,  for  case  2  in  Table  II.  Shown  are 
the  2D  primary  mode,  kz  =  0.2,  ky  =  0,  the  3D  primary  mode  kz  =  0.2,  ky  —  0.2,  and 
the  longest  wavelength  sheetwise  mode  kz  =  0,  ky  =  0.2.  Both  the  magnetic  and  velocity 
fields  evolve  in  essentially  the  same  way.  Rapid  initial  growth  followed  by  saturation  is 
seen  in  the  sheetwise  mode.  The  2D  and  3D  primary  modes  exhibit  a  burst  of  growth 
around  t  =  100,  followed  by  a  long  term  secular  decay.  The  crucial  feature  in  Figures 
11  and  12  is  the  rapid  growth  and  persistence  of  the  sheetwise  mode,  especially  in  the 
kinetic  energy.  These  modes  (i.e.,  all  those  with  kz  =  0)  eventually  dominate  the  system 
energetics.  Hence  we  see  a  burst  of  excitation  around  t  =  100,  followed  by  dominance  of 
the  purely  sheetwise  mode.  The  same  behavior  is  seen  in  higher  harmonics.  Preliminary 
calculations  at  higher  Reynolds  numbers  indicate  that  the  non-sheetwise  modes  do  not 
decay  as  quickly.  Consequently,  we  conjecture  that  the  relaminarization  is  a  dissipative 
effect  which  will  be  less  dominant  at  higher  Reynolds  numbers.  We  will  address  this  more 
fully  in  a  subsequent  paper. 

Some  morphological  features  of  the  evolution  of  the  secondary  instability  are  more 
easily  shown  in  configuration  space.  Figures  13  through  17  show,  in  the  (x,y,  z  =  0) 
plane,  the  time  evolution  of  contours  of  |B|2  for  run  1  in  Table  II.  Two  periods  in  x  are 
shown.  Figure  13  shows  these  contours  at  t  =  44,  slightly  before  the  rapid  development 
seen  in  figure  10.  In  this  figure,  the  leftmost  group  of  contours  shows  the  position  of  a 
magnetic  O-point,  the  next  set  a  magnetic  X-point,  and  so  forth.  A  sort  of  kinking  is 
evident  in  y,  with  the  O-points  somewhat  more  susceptible  than  the  X-points.  This  is  the 
basic  form  of  the  secondary  mode.  Figure  14  shows  the  increasingly  nonlinear  character 
of  the  disturbance.  Figure  15  shows  the  X-  and  O-points  beginning  to  interact.  Figures 
14  and  15  exhibit  the  continuing  formation  of  more  and  more  structure  in  the  sheetwise 
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direction.  Figures  16  and  17  show  later  stages  in  the  evolution  of  the  system.  It  of 
particular  interest  that  the  field  gradients  evolve  toward  being  in  the  sheetwise  direction, 
whereas  initially  they  were  in  the  fieldwise  direction.  Implicit  in  this  development  is  the 
breakdown  of  the  familiar  laminar  structures  associated  with  magnetic  reconnection,  e.g., 
the  magnetic  islands,  electric  current  filaments,  and  the  flow  pattern  shown  in  Figure  3. 
Little  topological  change  is  seen  once  the  contours  of  |B|2  adopt  the  pattern  shown  in 
Figure  17. 

Its  important  to  note  that  we  have  not  allowed  for  the  development  of  subharmonic 
modes,  *.e.,  the  perturbations  have  set  the  spatial  scales  in  the  fieldwise  and  sheetwise  di¬ 
rections.  In  random  noise  initiated  calculations  that  we  have  performed,  we  have  identified 
two  modes,  similar  to  modes  identified  for  velocity  free  shear  layers  41 .  (i.)  a  translational 
mode,  seen  in  figure  13,  and  (ii.)  a  pairing  mode,  which  occurs  when  subharmonic  in¬ 
teractions  are  possible.  This  last  mode  is  the  3D  analog  of  the  coalescence  mode,  and 
it  resembles  the  one  identified  in  low  Reynolds  number  experiments  38  38 .  The  presence 
of  these  two  modes  indicates  that,  as  in  hydrodynamical  problems,  the  route  by  which  a 
system  transitions  to  turbulence  i6  not  unique,  and  it  depends  heavily  on  initial  conditions. 

V.  CONCLUSIONS 

The  key  result  of  our  simulations  is  that  the  saturated  2D  reconnection  state  is  unsta¬ 
ble  to  3D  perturbations.  This  secondary  instability  is  non-dissipative  in  character,  and  in 
general  it  grows  at  a  rate  which  exceeds  the  reconnection  rate.  The  secondary  mode  growth 
rate  is  relatively  independent  of  the  dissipative  spreading  of  the  neutral  sheet,  at  least  on 
the  time  scales  over  which  our  simulations  were  performed.  When  the  secondary  mode 
becomes  nonlinear,  the  system  undergoes  a  turbulent  transition  to  a  three-dimensional 
quasi-steady  state.  During  the  late  nonlinear  phase,  the  sheetwise  modes  in  the  magnetic 
field  and  the  velocity  field  eventually  dominate  the  system. 

Hence,  it  is  legitimate  to  inquire  into  the  validity  of  the  standard  reconnection  sce¬ 
narios.  Our  work  suggests  that  the  conventional  current  sheet,  vortex  quadrupole  system 
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is  difficult  to  sustain  when  the  third  spatial  dimension  is  taken  into  account.  Therefore,  a 
steady  state  involving  such  a  configuration  is  unlikely.  The  observations  can  be  construed 
as  supporting  this  contention,  since  the  2D  laminar  model  consistently  underestimates  the 
energy  release  rates  required  to  account  for  phenomena  such  as  solar  flares.  Of  course,  we 
should  emphasize  that  so  far  we  have  investigated  only  neutral  sheets.  The  addition  of  a 
strong  field  component  in  the  sheetwise  direction  might  affect  the  evolution  of  an  electric 
current  sheet.  We  are  presently  addressing  this  issue. 
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Aj.  pendix  A:  Proof  of  Ideal  Stability  for  the  Primary  Equilibrium 

Consider  the  ideal  limit  of  equations  7,  i.e.,  let  R  — *  oo  and  Rm  — ►  oo  Furthermore, 
decompose  the  eigenvalue  as  u>  =  ac.  Then  equations  7  become 

0  =  —c{D2  -  (a2  +  p)w  +  {(D2B0)  -  B0[D2  -  (a2  +  02)}}bz  ,  (Ala) 

and 

cbz  =  —Bow,  (A16) 

where  we  have  set  Ma  =  1  and  Uo  =  0.  We  solve  equation  (Alb)  for  bz,  and  substitute 
the  result  in  equation  (Ala).  After  some  algebra,  this  reduces  to  the  form 

D[(c2  -  B2)Dw)  -  (a2  +  /?2)(c2  -  B2)w  =  0  .  (A2) 

Now  multiply  this  equation  by  the  complex  conjugate  of  w ,  w* ,  and  integrate  the  result  in 
the  crossfield  direction  using  free  slip  boundary  conditions.  After  some  manipulation,  we 
have  the  result: 

/+°V  -  Bl)[\Dw\*  +  (a2  +  02)M 2}dz  =  0.  (A3) 

J  —  OO 
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Finally,  separate  equation  (A3)  into  its  real  and  imaginary  parts. 


-  B’)[|Dw|!  +  (a1  +  fi)\w\2\iz  =  0, 


(Ala) 


[\Dw\2  +  (a2  +  02)\w\2)dz  =  0, 


(Mb) 


where  cr  +  *c<  =  c.  Here  equation  (A4a)  is  the  real  part  of  equation  (A3),  and  equation 
(A4b)  is  the  imaginary  part.  For  equation  (A4b)  to  be  true,  one  of  the  following  must  be 
true;  cr  =  0,  c<  =  0,  or  w(z)  =  0.  We  are  not  interested  in  the  final  option,  hence  one  of  the 


first  two  must  be  true.  An  inspection  of  equation  (A4a)  shows  that  the  equality  cannot  be 
true  if  cr  —  0,  and  hence  we  conclude  that 


Ci  =  0,  (.45) 

Since  u  =  ac  =  a(cr  +  *c,-),  and  the  primary  instabilities  vary  in  time  as  e~lut,  the  primary 
equilibrium  is  stable  in  the  ideal  limit. 


Appendix  B:  Numerical  Algorithm 

We  solve  equations  1  together  with  the  boundary  conditions  with  a  semi-implicit  spec¬ 
tral  scheme.  The  equations  are  discretized  in  the  two  periodic  directions,  x  and  y,  with 
a  Fourier  pseudospectral  scheme.  In  the  cross-sheet  direction,  z,  a  Fourier  collocation 
scheme  is  employed  together  with  a  hyperbolic  mesh  stretching.  The  method  extends  the 
algorithm  of  Cain  et  al.  42 ,  with  sine  and  cosine  functions  used  in  z,  together  with  a  cotan¬ 
gent  mapping  for  expansion  on  an  infinite  interval  (see  Canuto  et  al.  45 ,  pp  74-75).  The 
nonlinear  terms  are  advanced  in  time  with  a  third  order  Runge-Kutta  method,  while  the 
vertical  diffusion  is  advanced  with  a  Crank-Nicolson  scheme.  A  backward  Euler  pressure 
correction  is  performed  at  each  time  level.  The  same  type  of  correction  is  performed  at 
each  level  to  ensure  the  solenoidality  of  the  magnetic  field  44  45 ,  i.e.  we  define  a  scalar 
function  4>  such  that 

VV  +  V-B  =  0.  (Bl) 
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The  magnetic  field  is  corrected  after  the  advection-diffusion  step  so  that 


B'  =  B  +  V<£,  (B2) 

and  hence 

V  •  B'  =  0.  (£3) 


22 


References 


1.  E.  R.  Priest,  in  Solas  Flare  Magnetohydropdynanics ,  ed  E.  R.  Priest,  (Gordon  and 
Breach:  New  York,  1981),  p  139. 

2.  H.  R.  Strauss,  Phys.  Fluids  29,  3668  (1986). 

3.  H.  R.  Strauss,  Astrophys.  J.  326,  412  (1988). 

4.  D.  Biskamp  and  H.  Welter,  Phys.  Fluids  B  1,  1964  (1989). 

5.  W.  H.  Matthaeus  and  S.  L.  Lamkin,  Phys.  Fluids  29,  2513  (1986). 

6.  D.  Deeds  and  G.  Van  Hoven,  J.  Plasma  Phys.  42,  269  (1989). 

7.  Bayly,  B.  J.,  Orszag,  S.  A.,  and  Herbert,  Th.,  Ann.  Rev.  Fluid  Mech.  20,  359  (1988). 

8.  Herbert,  Th.,  Ann.  Rev.  Fluid  Mech.  20,  487  (1988). 

9.  Orszag,  S.  A.,  and  Patera,  A.  T.,  J.  Fluid  Mech.  128,  347  (1983). 

10.  Orszag,  S.  A.,  and  Kells,  L.  C.,  J.  Fluid  Mech.  96,  159  (1980). 

11.  Davey,  A.,  DiPrima,  R.  C.,  and  Stuart,  J.  T.,  J.  Fluid  Mech.  31,  17  (1968). 

12.  Kelly,  R.  E.,  J.  Fluid  Mech.  27,  657  (1967). 

13.  R.  W.  Metcalfe,  S.  A.  Orszag,  M.  E.  Brachet,  S.  Menon,  a^d  J.  J.  Riley,  J.  Fluid 
Mech.  184,  207  (1987). 

14.  Furth,  H,  Killeen,  J.  and  Rosenbluth,  M.  N.,  Phys.  Fluids  6,  459  (1963). 

15.  Rutherford,  P.  H.,  Phys.  Fluids  16,  19003  (1973). 

16.  White,  R.  B.,  Monticello,  D.  A.,  Rosenbluth,  M.  N.,  and  Waddell,  B.  V.,  Phys.  Fluids 
20,  800  (1977). 

17.  D.  Montgomery,  in  Spectral  Methods  for  Partial  Differential  Equations ,  eds  R.  G. 
Voight,  D.  Gottleib,  and  M.Y.  Hussaini  (Siam,  Philadelphia,  1984). 

18.  H.  B.  Squire,  Proc.  Roy.  Soc.  A  142,  621  (1933). 


23 


19.  D.  Montgomery,  Physica  Scripts  T2:2,  506  (1982). 

20.  J.  T.  Stuart,  Proc.  Roy.  Soc  A  221,  189  (1954). 

21.  E.  P.  Velikhov,  Sov.  Phys.  JETP  36,  848  (1959). 

22.  P.  R.  Nachtsheim  and  E.  Reshotko,  NASA  Tech.  Note  D-3144,  (1965). 

23.  Drazin,  P.  G.  and  Reid,  W.  H.,  Hydrodynamic  Stability  (Cambridge,  New  York, 
1984),  P156. 

24.  R.  B.  Dahlburg,  T.  A.  Zang,  D.  Montgomery,  and  M.  Y.  Hussaini,  Proc.  Nat.  Acad. 
Sci.  USA  80,  5798  (1983). 

25.  A.  Bondeson  and  J.  R.  Sobel,  Phys.  Fluids  27,  2028  (1984). 

26.  G.  Einaudi  and  F.  Rubini,  Phys.  Fluid  B  1,  2224  (1989). 

27.  R.  B.  Dahlburg,  T.  A.  Zang,  and  D.  C.  Montgomery,  J.  Fluid  Mech.  169,  71  (1986). 

28.  J.  T.  Stuart,  J.  Fluid  Mech.  4,  1  (1958). 

29.  R.  D.  Parker,  R.  L.  Dewar,  and  J.  I .  Johnson,  Phys.  Fluids  B  2,  509  (1990). 

30.  F.  Malara,  P.  Veltri,  and  V.  Carbone,  Phys.  Fluids  B  3,  1801  (1991). 

31.  R.  B.  Dahlburg  and  J.  M.  Picone,  Phys.  Fluids  B  1,  2153  (1989). 

32.  J.  M.  Picone  and  R.  B.  Dahlburg,  Phys.  Fluids  B  3,  29  (1991). 

33.  D.  C.  Montgomery,  in  Lecture  Notes  in  Turbulence,  eds.  J.  R.  Herring  and  J.  C. 
McWilliams  (World  Scientific,  Teaneck,  NJ,  1989),  p  75. 

34.  S.  A.  Orszag  and  Y.  -H.  Pao,  in  Turbulent  Diffusion  in  Environmental  Pollution,  ed. 
F.  N.  Frenkel,  and  R.  E.  Munn,  (Academic  Press,  New  York,  1974),  p  225. 

35.  D.  Schnack  and  J.  Killeen,  J.  Comp.  Phys.  35,  110  (1980). 

36.  S.  A.  Orszag,  J.  Fluid  Mech.  50,  689  (1971). 

37.  W.  H.  Matthaeus  and  D.  Montgomery,  J.  Plasma  Phys.  25,  11  (1981). 


24 


38.  W.  Gekelman  and  H.  Pfister,  Phys.  Fluids  31,  2017  (1988). 

39.  W.  Gekelman,  H.  Pfister,  and  J.  R.  Kan  J.  Geophys.  Res.  98,  3829  (1991). 

40.  I.  B.  Bernstein,  E.  A.  Frieman,  M.  D.  Kruskal,  and  R.  M.  Kulsrud,  Proc.  Roy.  Soc. 
(London)  A  244,  17  (1958). 

41.  R.  T.  Pierrehumbert  and  S.  E.  Widnall,  J.  Fluid  Mech.  114,  59  (1982). 

42.  Cain,  A.  B.,  Ferziger,  J.  H.,  and  Reynolds,  W.  C.,  J.  Comp.  Phys.  56,  272  (1984). 

43.  C.  Canuto,  M.  Y.  Hussaini,  A.  Quarteroni,  and  T.  A.  Zang,  Spectral  Methods  in  Fluid 
Dynamics  (Springer- Verlag,  New  York,  1988). 

44.  J.  U.  Brackbill  and  D.  C.  Barnes,  J.  Comp.  Phys.  35,  426  (1980). 

45.  Stephens,  A.  B.,  Bell,  J.  B.,  Solomon,  J.  M.,  and  Hackerman,  L.  B.,  J.  Comp.  Phys. 
53,  152  (1984). 


25 


TABLE  I. 


Case 

R  =  Rm 

0 

72  d 

lid 

7 

1 

0.2 

.042485 

.023348 

2 

0.2 

.043658 

.029808 

1 

3 

100 

0.2 

.033066 

.025587 

4 

100 

0.5 

.033066 

.007061 

l 

5 

100 

1.0 

.033066 

STABLE 

u 

6 

100 

1.5 

.033066 

STABLE 

7 

100 

2.0 

.033066 

STABLE 

.006 

8 

100 

2.5 

.033066 

STABLE 

-.023 

9 

200 

0.2 

.026636 

.020737 

.051 

10 

200 

0.5 

.026636 

.007288 

.062 

11 

200 

1.0 

.026636 

STABLE 

.059 

12 

200 

1.5 

.026636 

STABLE 

.051 

13 

200 

2.0 

.026636 

STABLE 

.039 

14 

200 

2.5 

.026636 

STABLE 

.023 

15 

200 

3.0 

.026636 

STABLE 

.004 

16 

270 

3.5 

.026636 

STABLE 

-.017 

17 

400 

0.2 

.020551 

.015915 

.051 

18 

1000 

0.2 

.013818 

.010528 

.048 

TABLE  II. 


Case 

a  =  0 

R  =  Rm 

62  d 

es  d 

72  d 

73  d 

1 

■ZH 

50 

.02 

.001 

jTS 

.02938 

2 

sm 

100 

.02 

.001 

■ 

.02559 

3 

0.2 

200 

.02 

.001 

.02074 
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Figure  1.  Primary  unstable  eigenmode  vs  z,  the  cross-sheet  direction.  All  modes  have 
a  =  .2,R  =  Rm  =  200.  a.  0  =  0,  2D  primary  unstable  mode;  b.  0  =  .2,  3D  primary 
unstable  mode;  c.  0  =  .5,  3D  primary  unstable  mode;  d.  0  =  2.,  a  nonpropagating 
3D  damped  mode. 


Figure  3.  Contour  plots  in  the  x  —  y  plane  exhibiting  the  structure  of  the  2D  saturated 
state  for  case  9  of  Table  I:  a.  magnetic  field,  b.  sheetwise  electric  current  density  and, 
c.  velocity  stream  function. 


ENERGIES 


TEARING  LAYER  ( ALPHA=BETA= . 2 , R=RM=200) 


Figure  4.  Time  evolution  of  the  logarithms  of  the  two  and  three-dimensional  total  energies 
for  case  9  in  Table  I.  Note  the  existence  of  a  well  defined  period  of  linear  growth  of 
the  3D  modes  after  the  2D  modes  saturate. 
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3D  GROWTH  COMPONENTS 


3D  TEARING  (ALPHfl=BETA=.2,R=RM=1000) 


Figure  8.  Time  evolution  of  Tij,T23,i?3  and  7  (3d  growth  rate)  for  case  18  of  Table  I. 
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Figure  9.  Time  evolution  of  6W  and  the  logarithms  of  the  2d  and  3d  total  energies  for 
case  9  of  Table  I. 
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Figure  10.  Time  evolution  of  the  enstrophies  for  cases  1,  2,  and  3  of  Table  II.  The  solid 
line  is  case  1,  the  dashed  line  is  case  2,  and  the  mixed  line  is  case  3. 
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VELOCITY  HARMONICS 
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Figure  11.  Time  evolution  of  velocity  modal  energies:  primary  modes  for  case  2  of  Table 
II. 
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Figure  12.  Time  evolution  of  magnetic  modal  energies:  primary  modes  for  case  2  of  Table 

II. 


38 


1050 


nnecTic  riEio  n«chituoc  oh  ncutrul  plane 


Two  periods  in  x  are  shown. 


f1A6NETlC  FIELD  ftH6N I TUDE  ON  NEUTRAL  PLANE 


>r 


► 


©  o  e  o  e  o 

a  o  c  o  o  o  o 

3  0  N  7  U)  O  Q 

OQOOOO*- 
>-0  000  00 

Z  *****  * 

0000000 

u 


40 


Figure  14.  Contour  plots  of  |B|2  at  t  =  133  in  the  ( x,y,z  =  0)  plane  for  case  1  of  Table 
II.  Two  periods  in  x  are  shown. 
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II.  Two  periods  in  x  are  shown. 
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II.  Two  periods  in  x  are  shown. 
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II.  Two  periods  in  x  are  shown. 


